Accounting for adsorption and desorption in Lattice Boltzmann simulations 
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We report a Lattice-Boltzmann scheme that accounts for adsorption and desorption in the cal- 
culation of mesoscale dynamical properties of tracers in media of arbitrary complexity. Lattice 
Boltzmann simulations made it possible to solve numerically the coupled Navier-Stokes equations 
of fluid dynamics and Nernst-Planck equations of electrokinetics in complex, heterogeneous media. 
Associated to the moment propagation scheme, it became possible to extract the effective diffusion 
and dispersion coefficients of tracers, or solutes, of any charge, e.g. in porous media. Neverthe- 
less, the dynamical properties of tracers depend on the tracer-surface affinity, which is not purely 
electrostatic, but also includes a species-specific contribution. In order to capture this important 
feature, we introduce specific adsorption and desorption processes in a Lattice-Boltzmann scheme 
through a modified moment propagation algorithm, in which tracers may adsorb and desorb from 
surfaces through kinetic reaction rates. The method is validated on exact results for pure diffusion 
and diffusion-advection in Poiseuille flows in a simple geometry. We finally illustrate the importance 
of taking such processes into account on the time-dependent diffusion coefficient in a more complex 
porous medium. 

PACS numbers: 47.00.00, 47.ll.-j, 47.55.dr, 47.56.+r, 47.61.-k, 47.70.Fw, 47.15. Rq, 47.80.-v 



I. INTRODUCTION 

The dynamical properties of fluids in heterogeneous 
materials offers a great challenge and have implications 
in many technological and environmental contexts. In- 
herently multi-scale in time and space, mesoscale prop- 
erties such as the diffusion or dispersion coefficient re- 
flect the nano-to-microscopic geometry of the media and 
the inter-atomic interactions between flowing particles, 
or tracers, and surface atoms. Experimentally, informa- 
tion about the microstructure of porous media can be ex- 
tracted from diffusion measurements by pulsed gradient 
spin echo nuclear magnetic resonance (PGSE-NMR) [1- 
3]. At short times, the dynamics of a pulse of tracers is 
connected to the geometry of the porous medium at the 
pore scale. At longer times, macroscale properties such 
as porosity and tortuosity come into play. Theoretically, 
stochastic approaches have been an important support 
to the understanding of the underlying phenomena [2-4] , 
and have been used recently to show that adsorption and 
desorption processes may strongly modify the short and 
long time dynamics of the tracers [5? ]. 

In numerous if not all practical situations involving 
particle diffusion and advection, the carrier fluid is in 
contact with confining walls where adsorption may oc- 
cur. These processes depend on the chemical nature of 
the solute, which explain why particles with the same 
charge can diffuse in the same medium with large differ- 
ences in their effective diffusion coefficients. This species- 



dependent affinity is at the heart of all chromatographic 
techniques used in analytic and separation chemistry [7] . 
It also plays a crucial role in the dissemination of toxic or 
radioactive pollutants in the environment, and conversely 
in remediation strategies. Recently, the great interest 
for nanofluidic devices and for the transport in hetero- 
geneous porous media has also raised the issue of the 
relevance of models which do not take into account these 
sorption processes. Moreover, it was recently shown that 
stochastic resonance between these processes and some 
external field may be of practical importance, e.g. for 
molecular sorting [5, 8]. 

At the mesoscale, the dynamics of particles in a fluid 
can be described by the continuity equation 



fitp(r,t) = -V-J(r,t), 



(1) 



where p(v, t) is the one-particle density at position r and 
time t, and J is the particle flux, which is a function of the 
velocity field of the carrier fluid, the bulk diffusion coeffi- 
cient Df, of the particles and, if any, their charge and the 
local electric field arising from their environment. If ad- 
sorption is taken into account, solid-liquid interfaces lo- 
cated at r have a surface concentration T (r, t) (length -2 ) 
that evolves with time according to: 



d t T (r,t) = -k d T (r,t)+k a p(r,t), 



(2) 
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where k a (length-time -1 ) and kd (time -1 ) are kinetic ad- 
sorption and desorption rates. We assume that the trac- 
ers (the solutes) neither diffuse into the solid phase (even 
though that process can easily be accounted for using our 
algorithm) nor dissolve the surfaces [9]. 
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The time-dependent diffusion coefficient D(t) and the 
dispersion coefficient K of the tracers can be investigated 
by following the spreading of a tracer pulse in the fluid. 
This would amount to solving Eqs. 1 and 2, e.g. with 
a finite element method, for all possible initial condi- 
tions, which is computationally intractable for complex 
systems such as heterogeneous porous media. An alter- 
native is to deduce K and D(t) from the tracer velocity 
auto-correlation function (VACF) following [10-13]: 

D 1 (t) = I Z 1 (t 1 ) dt', (3) 
Jo 

/•OO 

K 7 = / (Z 1 (t) - Z. ( (oo)) dt, (4) 
Jo 

where the VACF in the direction 7 e {x, y, z} is 

Z 7 (t) = (v 7 (0)v 7 (t)). (5) 

At long times Z 1 (00) = v 2 with v 1 the average velocity 
of the flow. The issue of averaging over initial conditions 
in Eq. 5 can be handled elegantly and efficiently using 
the moment propagation method [11, 14, 15], which was 
recently extended to charged tracers [13, 16, 17]. 

In order to compute the VACF of the tracers from 
Eq. 5, one has to keep track of their velocity. For this 
purpose, we use the underlying dynamics of the fluid 
given by Eq. 1, which does not rely on the velocity of 
individual particles, but on the one-particle solvent den- 
sity. Moreover, the simulation of heterogeneous multi- 
scale media requires a numerically efficient method. The 
lattice-Boltzmann (LB) method [18-24] offers a conve- 
nient framework to deal with such situations. In the 
LB approach, the fundamental quantity is a one-particle 
velocity distribution function fi (r, t) that describes the 
density of particles with velocity Cj, typically discretized 
over 19 values for three-dimensional LB, at a node r, ei- 
ther fluid or solid, of a lattice of spacing Ax, and at a 
time t discretized by steps of At. The dynamics of the 
fluid are governed by transition probabilities of a parti- 
cle moving in the fluid from one node to the neighboring 
ones: 

fi (r + c, At, t + At)= fi (r, t) + A t (r, t) , (6) 

where Aj, the so-called collision operator, is the change 
in fi due to collisions at lattice nodes. This LB equation 
recovers the fluid dynamics of a liquid and the moments 
of the distribution function are related to the relevant hy- 
drodynamic variables. The reader is referred to Refs. [22] 
and [24] for reviews of the method. 

II. ALGORITHM 

In order to compute the dynamical properties of trac- 
ers evolving in a fluid described by the LB method, 
i.e. to solve Eqs. 3-5, we use the moment propagation 
method [10, 11]. The moment propagation method has 



many advantages relative to other dispersion methods, 
the first of which is that it relies on the same ground as 
LB. It allows for the propagation of any moment of the 
particle velocity distribution function fi (r, t) [11, 14, 15], 
once a stationary flow state is reached. In the moment 
propagation algorithm, any quantity P (r, t) can be prop- 
agated between fluid nodes. This quantity will be mod- 
ified by adsorption and desorption processes. In their 
absence, P(r,t + At) = P* (r, t + At) with: 

P* (r, t + At) = i p ( r - c * A ^ *) Pi ( r - c * A ^ *)] 

i 

+P(r,()|l-^,(r,f)J, (7) 

where the first sum runs over all discrete velocities con- 
necting adjacent nodes. The probability of leaving node 
r along the direction Cj is noted pi (r, i). The last term 
in Eq. 7 represents the fraction of particles that did not 
move from r at the previous time-step. The expression 
for pi, which is central in the algorithm, depends on the 
nature of the tracers. It is given by: 

Pi(*,t) = -< tt -u>i + — Q, (8) 
p(r,t) 2 

where the first two terms account for advection and are 
obtained by coupling the tracer dynamics to that of the 
fluid evolving according to the LB scheme. The weigths 
u>i are constants depending upon the underlying LB lat- 
tice. The last term describes diffusive mass transfers. 
The dimensionless parameter A determines t he bulk dif- 
fusion coefficient D = \c 2 s At/ A, with c s = ^Jkj^Tjrn the 
sound velocity in the fluid. It also determines the mo- 
bility of tracers under the influence of chemical potential 
gradients (including the electrostatic contribution) which 
are accounted for in the Q term as described in Ref. [16]. 
For neutral tracers, Q = 1. 

We now introduce a new propagation scheme in order 
to account for adsorption and desorption at the solid- 
liquid interface. While Eq. 7 still holds for nodes r which 
are in the fluid but not at the interface, for the fluid 
interfacial nodes we define a new propagated quantity 
-Pads (r, t) associated with adsorbed particles: 

P ads (r, t + At) = P (r, t) Pa + P ads (r, t) (1 - p d ) , (9) 

where p a = k a At/ Ax is the probability for a tracer lying 
at an interface to adsorb, and pd — kdAt is the probabil- 
ity for an adsorbed, immobile tracer to desorb. There is 
no restriction in the definition of k a and kd so that they 
may depend on geometrical considerations and on the 
local tracer density. Finally, the evolution of the propa- 
gated quantity associated with free tracers now includes a 
term accounting for the desorption of adsorbed particles: 

P(r,t + At) = P*(r,t + At) + P ads (r,t)p d . (10) 
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where P* is still given by Eq. 7. In order to compute 
the VACF of the tracers, one propagates as P the prob- 
ability to arrive at position r at time t, weighted by the 
initial velocity of tracers. Thus one needs to initialize, 
for each direction 7, a propagated quantity according 
to the Maxwell-Boltzmann distribution. The Boltzmann 
weights for solid (S), fluid (J 7 ) and interfacial (I C J 7 ) 
nodes read: 



'0 



-/3A M " ds (r) 



for r e S, 

for re T \X, (11) 
/ Z for r e X, 



where e _ ^ AM ° *W = k a / (kdAx) corresponds to the sorp- 
tion free energy for interfacial tracers, /3 = l/fegT with 
ks the Boltzmann's constant and T the temperature, and 
Z is the partition function of the tracers. The excess 
chemical potential [f x (r) includes in the case of trac- 
ers with charge q a mean-field electrostatic contribution 
qip (r) with ip the local electrostatic potential. The VACF 
is then simply given as in the no sorption case [16] by: 



z 7 (t) = ^p(r,<) X nr 

r \ i 



Pi (r, t) c, 7 



(12) 



III. VALIDATION 

In order to validate this new scheme on exact results, 
we consider the diffusion and dispersion of tracers in a 
slit pore, i.e. between two walls at positions x — and 
x = L. The time-dependent diffusion coefficient D (t) in 
the direction normal to the wall is given by [? ]: 



D b 



1 



(2k a + k d L) 

2k d (k d 



x C 



k d L 

X 2 
s) sinh k 



X 3 ((kd + s) cosh k + kaX sinh k) 



(13) 



where s is the Laplace conjugate of time t, L 1 is the 
inverse Laplace transform, \ = ys/Db and K = X^/2. 

In Fig. 1, we compare the exact time-dependent dif- 
fusion coefficient D(i) of Eq. 13 with the one extracted 
from Lattice-Boltzmann simulations for different sorp- 
tion strength. This last quantity is defined by the frac- 
tion of adsorbed tracers. Unless otherwise stated, all 
simulations are performed within a slit pore of width 
L = 100 Ax and a bulk diffusion coefficient Db = 
10~ 2 Ax 2 /At. These values are chosen very conserva- 
tive, since the LB method is known to be efficient even 
in narrow slits (even for L < 10 Ax) and for a wide range 
of magnitudes in the diffusion coefficients [15, 21]. L and 
Db therefore account for a negligible part of the differ- 
ence with exact results. We can thus purposely assess 
the effect of the new algorithm only. Simulations in the 
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FIG. 1. (Color online) The time-dependent diffusion coeffi- 
cient, D(t), normalized by the bulk diffusion coefficient, Db, 
of neutral tracers in a slit pore, as extracted from Lattice- 
Boltzmann simulations using our new scheme (symbols) and 
from the reference exact solution of Eq. 13 (lines). Several 
fractions of adsorbed tracers, or sorption strength, are pre- 
sented: (black, circles) %; (red, squares) 1 %; (green, trian- 
gle up) 6 %; and (blue, triangle down) 60 %. 



presence of sorption reported in Fig. 1 correspond to a 
fixed sorption rate k a = 10 _1 Ax/ At and varying des- 
orption rates kd = 10 _1 , 10~ 2 and 10 _3 Af _1 resulting 
in the increasing sorbed fractions indicated on the figure. 
Excellent agreement is found, for all sorption strengths. 

At the initial time, D(t = 0) is given by the frac- 
tion of adsorbed tracers. At intermediate times, sorp- 
tion/desorption processes significantly decrease the slope 
of D(t) 7 as the partial immobilization at the surface slows 
down the exploration of the pore. We have shown re- 
cently that for this range of parameters, this slope is 
given by kd/(l + k a L/2Db) [? ]. In this illustration of 
diffusion between two walls, the confinement is total so 
that for sufficiently long times, the effective diffusion co- 
efficient decreases exponentially with time and tends to 
zero. 

In order to quantify the error on D(t) with respect to 
the exact result of Eq. 13, we plot in figure 2 the relative 
error on the slope and origin of the linear fit of log 10 D(t) 
for long times, i.e. for Dbt/L 2 > 0.5, as a function of k a 
and kd- In the whole range of k a and kd, the relative 
errors on the slope and origin remain under 5 % and 
2 %, respectively, which is highly satisfactory. 

The effect of a pressure gradient has also been studied 
on the same system. The resulting Poiseuille flow in- 
duces Taylor-Aris dispersion [25, 26] of the tracers with 
a dispersion coefficient K, which is known exactly in the 
presence of adsorption and desorption in the simple slit 
geometry [5]: 



K 
~D~b 



i + p: 



102y 2 + 18y + 1 D b 



2iJ 



(14) 



210(l + 2y) 3 L 2 k d (l + 2y) 
where P e = Lv/Db is the Peclet number and y = k a /kdL 
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FIG. 2. (Color online) Contour plots of the relative error of 
our Lattice Boltzmann scheme with respect to exact results, 
given by Eq. 13, on the slope (top) and origin (bottom) of the 
linear regression of the time-dependent diffusion coefficient of 
neutral tracers in a slit pore, as a function of the adsorption 
and desorption rates k a and kd. 



In Fig. 3, we compare the dispersion coefficient as cal- 
culated by LB with the exact results of Eq. 14 for vari- 
ous sorption strengths, as a function of the Peclet num- 
ber. Adsorption significantly increases the dispersion, as 
it slows down part of the tracers. The agreement between 
our scheme and exact results is excellent, even for strong 
adsorption and Peclet numbers above 100. 

As mentioned above, the time-dependence of the dif- 
fusion coefficient is a signature of the intrinsic geometric 
properties of a porous medium. The simplest model of 
such media consists of a compact face centered cubic (fee) 
lattice of spheres of radius R [27]. The porosity, i.e. the 
fraction of empty (or fluid) space, is 1 — tt/ (3\/2) ~ 26%. 
The unit cell contains 4 octahedral cavities of radius 
w 0.41i? connected by 8 smaller tetrahedral cavities of 
radius m 0.22R by a small channels of radius ~ 0.15i?. 
This fee lattice is illustrated in Fig. 4 for a lattice param- 
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FIG. 3. (Color online) Dispersion coefficient of neutral tracers 
in a slit pore in the direction of the flux, normalized by the 
bulk diffusion coefficient, as a function of the Peclet number, 
as extracted from our Lattice-Boltzmann scheme (symbols) 
and from the exact results (lines). Several fractions of ad- 
sorbed tracers, or sorption strength, are presented: (black, 
circles) %; (red, squares) 1 %; (green, triangle up) 6 %; and 
(blue, triangle down) 60 %. 
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FIG. 4. (Color online) Perspective view of a unit cell of a 
face-centered cubic packing of spheres, of lattice parameter 
100 Ax. Solid nodes are colored in blue. Interfacial fluid 
nodes, where adsorption processes may occur, are colored in 
red. Non-interfacial fluid nodes are in white. 



We report the time-dependent diffusion coefficient for 
this model porous medium in Fig. 5. At t = 0, the dif- 
fusion coefficient is again given by the fraction of free 
particles times their bulk diffusion coefficient. After a 
reduced time D^t/R 2 = 0.5, tracers have explored the 
whole porosity and the diffusion coefficient tends toward 
the effective diffusion coefficient. The time dependence 
is strongly influenced by adsorption/desorption, as in the 
slit pore case. It is thus essential to consider these phe- 
nomena when interpreting experimental measurement of 
effective and time-dependent diffusion coefficients. 
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FIG. 5. (Color online) Diffusion coefficient D(t) of a neu- 
tral tracer, normalized by the bulk diffusion coefficient Db in 
the face-centered cubic packing of spheres with radius R il- 
lustrated in Fig. 4, as a function of the reduced time. Several 
fractions of adsorbed tracers, or sorption strength, are pre- 
sented: (black, circles) %, i.e. without adsorption; (red, 
squares) 30 %; (green, triangle up) 79 %; and (blue, trian- 
gle down) 97 %. 



IV. CONCLUSION 

In summary, we have proposed a new scheme that 
accounts for adsorption and desorption in a generic 
Lattice-Boltzmann scheme, allowing for the calculation 
of mesoscale dynamical properties of tracers in media of 
arbitrary complexity. These processes are modelled by 
kinetic rates of adsorption and desorption taking place 
at interfacial, fluid nodes. The algorithm has been val- 



idated over a wide range of adsorption and desorption 
rates and Pcclct numbers in the slit pore geometry where 
exact results are available. Finally, we have shown on a 
more complex porous medium that adsorption and des- 
orption processes may not be neglected, as they strongly 
modify the short, intermediate and long time behaviors 
of the diffusion coefficient as well as the dispersion coef- 
ficient. In turn, this demonstrates that neglecting inter- 
actions with the surface in the interpretation of D(t) as 
a probe of the geometry of the porous medium, as mea- 
sured experimentally e.g. by PSGE-NMR, may lead to 
incorrect conclusions. This scheme may now be used in 
two ways. First, one could predict the effective diffusion 
coefficient in complex heterogeneous media for species 
with known adsorption and desorption rates (from ex- 
periments or molecular simulations). Conversely, from 
reference measurements of the time-dependent diffusion 
coefficient in controlled geometries, one could extract the 
adsorption and desorption rates k a and k^. Moreover, in 
the case of diffusion in the solid stationary phase, the 
method would allow us to relate the relevant diffusion 
constant to the shape of the elution profile. 
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